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An  investigation  to  improve  trajectory  prediction  using  Lagrangian  data  is  presented.  The  velocity  field  of 
a  data  assimilating  model,  EAS-16,  is  corrected  using  drifter  observations  taken  during  an  experiment  off 
Taiwan.  The  results  are  tested  using  another  independent  Lagrangian  data  set  provided  by  sonobuoys 
launched  in  the  same  area.  The  latter  have  instrument  chains  that  extend  well  into  the  water  column. 
Consequently  the  corrected  model  velocities  were  projected  into  the  water  column  in  order  to  calculate 
sonobuoy  trajectories.  The  drifter  and  sonobuoy  trajectories  both  show  two  distinct  regimes  in  the  con¬ 
sidered  area  of  approximately  1/2°  square.  One  regime  is  dominated  by  shelf  dynamics,  the  other  by 
meandering  of  the  Kuroshio,  with  a  sharp  boundary  dividing  the  two.  These  two  regimes  are  not  repro¬ 
duced  by  the  trajectories  of  the  EAS-16  model.  When  the  drifter  data  are  blended  with  the  model  veloc¬ 
ities,  synthetic  sonobuoy  trajectories  track  the  observed  ones  much  better,  and  the  two  regimes  are 
clearly  depicted.  Two  different  methods  for  the  velocity  reconstruction  are  tested.  One  is  based  on  a  var¬ 
iational  approach  and  the  other  on  a  normal  mode  decomposition.  Both  methods  show  qualitatively  sim¬ 
ilar  improvements  in  the  prediction  of  sonobuoys  trajectories,  with  a  quantitative  improvement  in  the 
total  rms  error  of  approximately  50%  and  25%,  respectively. 

©  2010  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Operational  forecasting  systems  are  based  on  two  complemen¬ 
tary  components,  a  monitoring  system  and  an  ocean  modeling  sys¬ 
tem.  The  model  ocean  state  is  routinely  corrected  using  the  data 
from  the  monitoring  system,  and  forecasts  are  provided  running 
the  model  forward  starting  from  a  corrected  (constrained)  initial 
state.  Forecast  skills  have  dramatically  increased  in  the  last  years, 
but  their  main  limiting  factor  may  well  be  related  to  the  density 
and  quality  of  the  observations  that  are  used  to  constrain  the  anal¬ 
ysis  and  re-initialize  the  model. 

For  many  applications,  there  is  now  a  significant  demand  for 
Lagrangian  products  from  operational  models.  Examples  are  search 
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and  rescue,  drifting  sensor  arrays,  and  mitigation  in  case  of  pollu¬ 
tants,  such  as  oil  spills.  Lagrangian  predictability,  i.e.  prediction  of 
particle  motion,  is  especially  demanding  for  a  number  of  reasons. 
Particle  trajectories  are  obtained  by  integrating  the  velocity,  so  that 
even  small  errors  in  the  forecasts  of  Eulerian  velocity  tend  to  accu¬ 
mulate  and  grow  (Griffa  et  al.,  2004).  Also,  particle  motion  is  often 
inherently  chaotic,  namely  it  exhibits  a  high  dependence  on  initial 
conditions  even  in  very  simple  flows  (Aref,  1984).  Thus,  even  a 
slight  difference  in  initial  conditions  in  space  and  time  can  result 
in  significantly  different  behaviors.  A  natural  avenue  to  improve 
trajectory  prediction  appears  to  be  tbe  assimilation  of  Lagrangian 
data  that  provide  direct  trajectory  information. 

In  the  last  decade,  a  number  of  schemes  for  using  Lagrangian  data 
for  forecasting  have  been  proposed.  Lagrangian  instruments  are 
floating  devices  (acting  as  proxies  for  fluid  particles)  that  provide 
information  on  their  positions  and  possibly  on  other  environmental 
parameters  at  discrete  time  intervals.  Velocity  information  can  be 
obtained  from  consecutive  position  observations,  provided  that 
the  time  interval  is  smaller  than  the  typical  Lagrangian  time  scale 
Ti,  namely  the  time  over  which  particle  velocity  is  self-correlated. 


Report  Documentation  Page 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 


1.  REPORT  DATE 

NOV  2010 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

Enhanced  estimation  of  sonobuoy  trajectories  by  velocity  reconstruction 
with  near-surface  drifters 

6.  AUTHOR(S) 


7.  PEREORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Research  Laboratory, Stennis  Space  Center, MS, 39529 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2010  to  00-00-2010 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

An  investigation  to  improve  trajectory  prediction  using  Lagrangian  data  is  presented.  The  velocity  field  of 
a  data  assimilating  model,  EAS-16,  is  corrected  using  drifter  observations  taken  during  an  experiment  off 
Taiwan.  The  results  are  tested  using  another  independent  Lagrangian  data  set  provided  by  sonobuoys 
launched  in  the  same  area.  The  latter  have  instrument  chains  that  extend  well  into  the  water  column. 
Consequently  the  corrected  model  velocities  were  projected  into  the  water  column  in  order  to  calculate 
sonobuoy  trajectories.  The  drifter  and  sonobuoy  trajectories  both  show  two  distinct  regimes  in  the 
considered  area  of  approximately  1/2  square.  One  regime  is  dominated  by  shelf  dynamics,  the  other  by 
meandering  of  the  Kuroshio,  with  a  sharp  boundary  dividing  the  two.  These  two  regimes  are  not 
reproduced  by  the  trajectories  of  the  EAS-16  model.  When  the  drifter  data  are  blended  with  the  model 
velocities  synthetic  sonobuoy  trajectories  track  the  observed  ones  much  better,  and  the  two  regimes  are 
clearly  depicted.  Two  different  methods  for  the  velocity  reconstruction  are  tested.  One  is  based  on  a 
variational  approach  and  the  other  on  a  normal  mode  decomposition.  Both  methods  show  qualitatively 
similar  improvements  in  the  prediction  of  sonobuoys  trajectories,  with  a  quantitative  improvement  in  the 
total  rms  error  of  approximately  50%  and  25%,  respectively. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OE 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

19 

Standard  Form  298  (Rev.  8-98} 

Prescribed  by  ANSI  Std  Z39-18 


180 


Y.  Chang  et  all  Ocean  Modelling  36  (2011)  179-197 


For  the  surface  ocean,  typically  varies  in  the  range  of  1-5  days 
(Bauer  et  al.,  2002;  LaCasce,  2008).  Various  methods  have  been 
suggested  for  both  velocity  reconstruction  (Toner  et  al.,  2001; 
Taillandier  et  al.,  2006a)  and  actual  assimilation  (Molcard  et  al., 
2003;  Taillandier  et  al.,  2006b;  Kuznetsov  et  al.,  2003).  The  former 
consists  of  improving  the  velocity  field  from  off-line  circulation 
model  output  by  blending  it  with  the  Lagrangian  data,  while  assim¬ 
ilation  implies  that  the  circulation  model  itself  is  corrected  (often  by 
sequential  re-initialization)  using  the  corrected  velocity  fields  and 
adequately  balanced  mass  fields.  Reconstruction  can  be  seen  as  a 
first  step  toward  assimilation,  but  it  is  also  valuable  per  se,  especially 
for  operational  purposes,  since  it  is  extremely  flexible  and  can  be 
applied  to  the  output  of  any  model  used  in  case  of  an  accident  or 
emergency,  even  when  the  model  does  not  have  a  full  assimilation 
capability. 

Two  main  approaches  have  been  followed  for  the  reconstruc¬ 
tion  and  assimilation  of  Lagrangian  data.  The  first  approach  is 
based  on  estimating  velocities  along  trajectories  as  the  ratio  be¬ 
tween  observed  position  differences  and  time  increments  (e.g., 
Hernandez  et  al.,  1995)  and  directly  using  these  velocities  to  cor¬ 
rect  the  model  results.  The  second  approach  introduces  an  obser¬ 
vational  operator  based  on  the  particle  advection  equation  and 
corrects  the  Eulerian  velocity  field  by  requiring  the  minimization 
of  the  difference  between  observed  and  modeled  trajectories 
(e.g.,  Molcard  et  al.,  2003).  These  approaches  have  been  imple¬ 
mented  using  several  methodologies,  including  optimal  interpola¬ 
tion  (01;  Molcard  et  al.,  2003;  Ozgbkmen  et  al.,  2003;  Molcard 
et  al.,  2005),  mode  decomposition  techniques  (Toner  et  al.,  2001), 
Kalman  filtering  (Ide  et  al.,  2002;  Kuznetsov  et  al.,  2003;  Salman 
et  al.,  2006),  variational  methods  (Kamachi  and  O’Brien,  1995; 
Taillandier  et  al.,  2006a;  Taillandier  and  Griffa,  2006;  Nodet, 
2006),  and  particle  filter  methods  (Salman,  2008b,a;  Krause  and 
Restrepo,  2009).  Apte  et  al.  (2008)  have  developed  a  Markov-chain 
Monte  Carlo  sampling  strategy  for  Lagrangian  data  assimilation 
that  eliminates  problems  with  the  Kalman  filter  approach  near 
saddle  points  in  the  flow  field.  Such  Bayesian  methods  should 
play  an  increasingly  important  role  in  assimilating  and  blending 
Lagrangian  and  Eulerian  data. 

While  these  methods  have  been  thoroughly  tested  using  syn¬ 
thetic  data  in  the  framework  of  numerical  models,  application  to 
in  situ  data  and  actual  operational  testing  are  just  beginning.  The 
variational  method  of  Taillandier  et  al.  (2006a)  has  been  applied 
to  Argo  float  data  assimilation  in  the  northwestern  Mediterranean 
Sea  (Taillandier  et  al.,  2006b,  2010)  and  to  drifter  data  reconstruc¬ 
tion  in  the  Adriatic  Sea  (Taillandier  et  al.,  2008).  An  application  of 
the  mode-decomposition  method  to  50  m  drogued  drifters  and  a 
primitive  equation  model  of  the  Gulf  of  Mexico  has  been  reported 
in  Toner  et  al.  (2001).  These  studies  show  significant  and  consis¬ 
tent  changes  in  the  ocean  circulation  and  in  the  Lagrangian  path¬ 
ways,  but  a  full  and  quantitative  evaluation  of  the  approaches 
using  independent  data  is  still  lacking.  This  study  provides  the  first 
such  assessment. 

Here,  we  address  this  issue  with  a  unique  data  set  collected 
from  an  exercise  off  Taiwan  in  October  2007.  The  data  are  com¬ 
posed  of  30  SVP  (Surface  Velocity  Program)  drifters  and  28  sonobu- 
oys  with  instrumented  chains  deployed  in  a  small  grid 
(approximately  1  /2“  square)  over  three  days  during  a  Littoral  War¬ 
fare  Advanced  Development  experiment  (LWAD07).  The  SVP  drift¬ 
ers  are  carefully  designed  to  follow  the  flow  field  at  approximately 
15  m,  while  the  sonobuoys  respond  to  the  flow  along  the  entire 
length  of  the  instrument  chain.  In  addition  to  this  data,  hindcasts 
for  the  experiment  are  provided  from  a  data  assimilating  model, 
the  Naval  Research  Laboratory  East  Asian  Seas  1/16°  ocean  model 
(EAS-16).  In  the  present  investigation,  the  velocity  field  is  recon¬ 
structed  using  two  different  methods  that  blend  the  output  of 
EAS-16  with  the  data  from  the  SVP  drifters  and  statistically  project 


the  velocity  correction  over  the  water  column.  The  corrected  veloc¬ 
ity  fields  are  then  independently  tested  using  the  sonobuoy  data. 
Synthetic  sonobuoy  trajectories  are  computed  using  an  appropri¬ 
ate  drag  model  applied  to  the  corrected  velocities,  and  they  are 
quantitatively  compared  with  the  observed  sonobuoy  trajectories. 
The  dense  array  of  Lagrangian  data  from  drifting  sensors  that  re¬ 
spond  to  currents  at  different  depths  and  the  near-operational  high 
resolution  data  assimilating  model  provide  an  opportunity  to  eval¬ 
uate  Lagrangian  predictability. 

Our  investigation  has  a  number  of  unique  and  novel  aspects.  It 
provides  an  example  of  a  truly  operational  application,  targeted  to 
predict  the  motion  of  sonobuoys  within  the  framework  of  an 
LWAD  exercise.  The  region  of  application  is  interesting  and  chal¬ 
lenging,  since  it  is  located  along  the  shelf  break  at  the  boundary 
of  the  Kuroshio,  encompassing  two  distinct  flow  regimes  in  a  rela¬ 
tively  small  area,  divided  by  a  geostrophic  front.  Predicting  the  ex¬ 
act  location  of  the  front  at  such  scales  is  challenging  for  a 
numerical  model.  Also,  the  use  of  two  different  methods  of  recon¬ 
struction  provides  an  interesting  opportunity  to  evaluate  the  util¬ 
ity  of  Lagrangian  data  blending  with  different  approaches,  finally 
and  most  importantly,  the  present  application  provides  a  first 
example  of  independent  testing  of  the  results,  since  the  reconstruc¬ 
tion  is  based  on  drifter  data  only,  while  the  testing  is  performed  using 
the  sonobuoy  trajectories.  Sonobuoys  not  only  respond  differently  to 
currents  than  drifters  do,  but  they  also  have  been  launched  at 
slightly  different  times  (order  of  one  day  difference)  from  the  drift¬ 
ers,  which  is  significant  for  Lagrangian  applications  characterized 
by  time  scales  of  1-2  days.  Positive  results  in  the  case  of  sonobuoy 
trajectory  prediction  are  expected  to  be  meaningful  also  for  other 
applications  that  involve  forecasts  of  floating  quantities  in  the 
upper  ocean.  Potential  applications  include  prediction  of  the  path¬ 
ways  of  pollutants  or  invasive  species  such  as  jelly  fish. 

The  paper  is  organized  as  follows.  The  LWAD07  experiment  and 
the  data  sources  are  described  in  Section  2.  In  Section  3,  the  two 
methods  used  to  blend  the  Lagrangian  data  with  EAS-16  are  dis¬ 
cussed  and  their  results  are  shown.  (More  details  on  the  methods 
are  provided  in  the  Appendices.)  The  assessment  of  the  recon¬ 
structed  velocities  using  sonobuoy  data  is  given  in  Section  4.  We 
conclude  with  a  summary  of  our  findings  and  a  discussion  of  the 
role  of  Lagrangian  data  in  enhancing  Lagrangian  forecasts. 

2.  Data  sources 

2.1.  LWAD07:  SVP  drifters 

The  data  set  we  are  using  for  the  reconstruction  consists  of  30 
SVP  drifters,  whose  configuration  followed  the  standards  of  the 
Global  Drifter  Program.  The  drifter’s  drogue  is  centered  at  15  m  be¬ 
neath  the  surface,  thus  ensuring  that  the  drifters  capture  near-sur- 
face  velocities.  Determining  drogue  loss  remains  challenging,  even 
when  submergence  measurements  are  provided,  as  in  this  case. 
Using  the  standard  procedure  of  equating  a  sharp  drop  in  the  sub¬ 
mergence  indicator  with  drogue  loss,  only  one  of  the  drifters  lost 
its  drogue  during  the  analysis  period.  We  have  chosen  to  include 
its  full  data  record  nonetheless,  since  its  statistics  and  behavior 
are  indistinguishable  from  the  others’. 

The  drifters  were  launched  during  the  window  of  October  8-12, 
2007  UTC  in  the  East  China  Sea.  The  deployment  area  straddles  the 
edge  of  the  continental  shelf  and  is  thus  near  the  Kuroshio,  a  strong 
western  boundary  current.  The  shelf  dynamics  are  dominated  by 
the  tides.  An  overview  of  the  geographical  and  topographical  con¬ 
text  of  the  launch  area  is  given  in  Fig.  1 ,  along  with  a  sample  model 
velocity  field.  All  but  one  of  the  drifters  reported  continuously  for 
at  least  the  first  two  weeks,  and  most  of  them  up  to  25  days,  the 
exception  being  a  drifter  that  died  after  just  one  day  in  the  water. 
GPS  positions  and  sea  surface  temperatures  were  recorded  at 
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Fig.  1.  Overview  of  drifter  launch  locations.  Purple  disks  indicate  launches  performed  en  route  to  the  standard  stations,  while  red  disks  indicate  the  four  standard  stations. 
Bottom  topography  is  displayed  in  the  color  field  to  illustrate  the  shelf  edge.  The  velocity  field  is  a  typical  EAS-16  model  field.  (For  interpretation  of  the  references  to  colour  in 
this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


approximately  fifteen  minute  intervals.  In  a  preprocessing  step, 
duplicate  entries  and  those  implying  unrealistically  large  velocities 
or  accelerations  were  deleted  from  the  records. 

Details  on  the  launch  times  and  on  time  coverage  during  the 
first  few  days  of  the  experiment,  that  are  the  main  focus  of  the 
present  application,  are  shown  in  Fig.  2(a).  The  first  six  drifters 
were  launched  around  22:00  October  8,  while  the  ship  was  en 
route  (purple  dots  in  Fig.  1 ),  while  the  other  drifters  were  deployed 
in  six  successive  deployments  of  four  instruments  each,  reusing 
four  standard  launch  locations  each  time  (red  dots  in  Fig.  1).  The 
six  successive  deployments  roughly  correspond  to  22:00  October 
9,  4:00  October  10,  15:00  October  10,  4:00  October  11,  16:00 
October  11,  and  4:00  October  12.  Each  event  of  four  launches  spans 
approximately  1.5  h,  except  for  one  event  which  includes  a  launch 
4  h  prior  to  the  rest.  Also,  one  event  had  two  southwest  launches 
but  none  from  the  northwest  corner. 

The  drifter  trajectories  corresponding  to  the  first  two  days  of 
the  experiment  (October  8-10,  2007)  are  shown  in  Fig.  3(a).  The 
motion  of  the  drifters  can  be  categorized  by  the  region  of  the  test 
site  in  which  they  were  deployed,  and  it  shows  the  presence  of 
two  very  distinct  regimes.  Drifters  deployed  in  the  northwestern 


(a)  Drifters  (b)  Sonobuoys 


Fig.  2.  Time  spans  of  (a)  drifters  and  (b)  sonobuoys  over  the  period  October  8-14, 
2007. 


182 


Y.  Chang  et  al./ Ocean  Modelling  36  (2011)  179-197 


(a)  In-Situ  drifters  (b)  In-Situ  Sonobuoys  (c)  EAS-16  drifters 


Longitude  Longitude  Longitude 


Fig.  3.  (a)  Observed  drifter  trajectories,  (b)  observed  sonobuoy  trajectories,  and  (c)  modeled  drifter  trajectories  using  EAS-16.  during  the  first  two  days  of  the  experiment. 
October  8-10  2007.  Colors  depict  trajectories  of  corresponding  drifters  in  panels  (a)  and  (c).  Black  circles  mark  initial  locations. 


corner,  characterized  by  shallower  water  and  shelf  dynamics,  show 
a  looping  pattern  with  almost  no  net  drift.  The  semi-diurnal  period 
of  the  looping  is  indicative  of  tidal  forcing.  The  drift  of  instruments 
deployed  in  the  southeastern  corner  of  the  test  site  is  influenced  by 
the  Kuroshlo  as  characterized  by  rapid  advectlon  to  the  northeast 
with  little  tidal  fluctuation.  Drifters  deployed  In  the  transition  zone 
between  the  two  corners  exhibit  characteristics  of  both  the  south¬ 
eastern  and  northwestern  corner  deployments.  These  drifters  gen¬ 
erally  start  slowly  looping  to  the  east  or  northeast.  Some  of  them 
are  pulled  closer  Into  the  influence  of  the  Kuroshio,  and  then  even¬ 
tually  they  get  swept  up  into  the  Kuroshio  and  move  rapidly  to  the 
northeast. 

2.2.  LWAD07:  CDMR  sonobuoys 

In  addition  to  the  SVP  drifters,  28  Compact  Deployable  Multi- 
static  Receiver  (CDMR)  sonobuoys  were  deployed  during  the 
LWAD07  experiment  from  October  8-11,  2007  in  the  same  target 
area  as  the  SVP  drifters.  The  CDMR  are  experimental  Navy  sonobu¬ 
oys  similar  in  physical  characteristics  to  the  AN/SSQ-101  ADAR 
sonobuoy  (Fig.  4).  They  consist  of  an  inflatable  toroid  shaped  sur¬ 
face  float  with  an  approximate  diameter  of  0.45  m,  a  length  of 
cable  and  bungee,  two  cylindrical  cable  packs,  and  a  large  horizon¬ 
tal  planar  hydrophone  array  (diameter  about  6  m)  set  to  a 
maximum  depth  of  27.4  m.  Each  sonobuoy  communicated  its  posi¬ 
tion  to  a  logging  station  via  an  Iridium  satellite  link. 


The  launch  times  and  life  spans  of  the  sonobuoys  are  summa¬ 
rized  in  Fig.  2(b).  The  average  operating  life  span  of  the  deployed 
CDMR  during  this  exercise  was  43  h  with  a  maximum  of  82  h. 
The  first  deployment  of  eight  sonobuoys  occurred  around  14:00 
October  8,  i.e.  approximately  7  h  before  the  first  drifter  deployment. 
Subsequent  deployments  span  the  period  from  22:00  October  8  to 
0:00  October  11.  Some  of  the  sonobuoys  were  launched  simulta¬ 
neously  with  drifter  deployments,  but  in  many  cases  there  is  no 
direct  correspondence  between  sonobuoy  and  drifter  launch  times. 
Moreover,  launch  locations  are  slightly  different. 

Sonobuoy  trajectories  are  shown  in  Fig.  3(b)  during  the  same 
time  interval  as  the  drifters  in  Fig.  3(a).  The  sonobuoy  trajectories 
show  a  similar  general  pattern  as  the  drifters,  characterized  by  the 
two  flow  regimes,  but  they  also  show  significant  differences. 
Sonobuoys  launched  in  the  northwestern  area  show  a  more  pro¬ 
nounced  southward  drift  than  the  drifters,  and  those  entrained  in 
the  Kuroshio  experience  a  somewhat  slower  drift.  The  differences 
are  likely  due  to  two  main  factors:  (a)  the  different  response  of 
sonobuoys  to  currents  at  different  depths  and  (b)  the  differences 
in  launch  times  and  locations.  Given  the  high  variability  of  the  flow 
and  the  inherent  chaotic  nature  of  Lagrangian  motion,  even  slight 
differences  in  the  initial  conditions  can  lead  to  significantly  differ¬ 
ent  trajectories. 

2.3.  The  EAS- 1 6  system 

The  Naval  Research  Laboratory  East  Asian  Sea  1/16°  ocean  mod¬ 
el  (EAS-16)  is  exercised  daily  using  atmospheric  forcing  provided 
by  the  Navy  Operational  Global  Atmospheric  Prediction  System 
(NOGAPS).  It  is  an  adaptation  of  the  operational  Navy  Coastal 
Ocean  Model  (NCOM;  Martin,  2000)  and  takes  its  boundary  forcing 
from  a  global  1/8°  NCOM  implementation.  These  are  primitive 
equation  models,  with  hydrostatic,  Boussinesq,  and  3D-incom- 
pressibility  approximations.  41  mixed-coordinate  layers  are  used, 
with  a  coordinates  from  the  surface  down  to  the  shelf  break 
around  137  m  depth  and  z  levels  below.  For  our  analysis,  the  out¬ 
put  has  been  interpolated  to  constant  depths  at  1,  2,  5,  10,  20,  40, 
60,  80, 100,  and  150  m.  Eight  tidal  components  (Kl,  01,  PI,  Ql,  K2, 
M2,  N2,  and  S2)  are  incorporated  through  tidal  potential  forcing 
from  NCOM.  The  model  assimilates  sea  surface  height  anomaly 
and  temperature  through  the  Modular  Ocean  Data  Assimilation 
System  (MODAS;  Fox  et  al.,  2002).  For  a  more  detailed  description 
of  the  model  and  a  discussion  of  the  Eulerian  validations  carried 
out,  see  Riedlinger  et  al.  (2006). 

While  EAS-16  is  employed  here  mainly  because  it  was  the  oper¬ 
ational  model  accompanying  the  LWAD07  experiment,  it  needs  to 
be  stated  that  EAS-16  is  not  what  the  Navy  considers  to  be  the 


y.  Chang  et  al.  /  Ocean  Modelling  36  (2011)  179-197 


183 


State-of-the-art  forecasting  system  anymore.  Models  with  1  km 
and  3  km  horizontal  resolutions  and  other  data  assimilations  sys¬ 
tems  have  been  used  in  the  LWAD  experiment  conducted  in 
2009.  Results  from  this  more  recent  application  will  be  reported 
in  a  future  study. 

Examples  of  drifter  trajectories  computed  using  the  EAS-16  out¬ 
put  are  shown  in  Fig.  3(c)  for  the  same  period  as  for  the  in  situ 
drifters  (Fig.  3(a)).  The  simulated  trajectories  are  initialized  at  the 
launch  locations  of  the  observed  drifters  and  are  computed  using 
a  4tb  order  Runge-Kutta  method.  Compared  with  Fig.  3(a),  the 
EAS-16  trajectories  show  a  more  homogeneous  drift  toward  the 
northeast  and  a  more  pronounced  looping  throughout  the  area, 
without  indicating  any  significant  difference  between  the  two  re¬ 
gimes,  on  the  shelf  and  within  the  Kuroshio.  This  suggests  that 
the  details  of  the  Kuroshio  meandering  are  not  correctly  repro¬ 
duced  in  the  model,  which  is  not  surprising  given  the  high  nonlin¬ 
earity  of  the  region  and  the  small  size  of  the  target  domain  relative 
to  the  resolution  of  the  model  and  of  the  assimilation  observing 
system. 

3.  Blending  Lagrangian  drifters  and  EAS-16  output 

In  order  to  improve  upon  the  Lagrangian  prediction  results  of 
EAS-16  cited  above,  we  employ  a  tool  that  has  proven  useful  for 
extending  Eulerian  predictability  horizons,  namely  data  assimila¬ 
tion.  Since  we  are  primarily  interested  in  Lagrangian  forecasts,  it 
is  natural  to  consider  Lagrangian  observations  and  how  best  to 
blend  them  with  the  model  output.  As  discussed  in  the  Introduc¬ 
tion,  correcting  off-line  velocities  is  different  from  performing  full 
assimilation  including  balancing  of  the  mass  variables  and  re-ini- 
tialization  of  the  model.  On  the  other  hand,  the  procedure  can  be 
seen  as  a  first  step  toward  assimilation,  and  it  has  the  advantage 
of  being  flexible,  portable,  and  easy  to  implement. 

Tbe  reconstruction  of  the  velocity  field  using  drifter  data  and 
EAS-16  output  is  performed  using  two  methods  based  on  different 
approaches.  One  method,  Lagrangian  variational  analysis  (LAVA),  is 
based  on  the  variational  approach  put  forward  by  Taillandier  et  al. 
(2006a),  where  model  velocities  are  corrected  by  minimizing  tbe 
distance  between  observed  and  modeled  drifter  positions.  Tbe 
other  method.  Normal  Mode  Analysis  (NMA),  is  based  on  estimat¬ 
ing  velocities  along  trajectories  and  optimally  blending  them  with 
model  velocities  (Toner  et  al.,  2001).  For  both  methods,  a  two  step 
strategy  is  adopted.  First,  the  drifter  data  are  used  to  correct  the 
EAS-16  velocity  at  a  10-m  depth,  approximately  corresponding  to 
the  drogue  nominal  depth  of  1 5  m.  The  velocity  correction  is  then 
projected  into  the  water  column  by  using  covariance  methods 
developed  from  the  original  unblended  EAS-16  output.  Starting 
with  the  20-m  model  level  instead  makes  little  difference,  as  there 
is  little  vertical  shear  at  these  depths  in  the  model.  Without  data 
blending,  the  rms  separation  distances  between  observed  and 
model  trajectories  at  10  m  and  at  20  m  after  24  h  differ  by  less  than 
0.05%  or  about  7  m. 

3.1.  Lagrangian  variational  analysis  (LAVA) 

Tbe  LAVA  approach  (Taillandier  et  al.,  2006a)  has  been  previ¬ 
ously  applied  to  Argo  float  data  assimilation  in  the  northwestern 
Mediterranean  Sea  (Taillandier  et  al.,  2006b,  2010)  and  to  drifter 
data  reconstruction  in  tbe  Adriatic  Sea  (Taillandier  et  al.,  2008). 
Some  details  of  the  method  are  given  in  Appendix  A,  while  the  spe¬ 
cific  configuration  used  here  is  discussed  below. 

Tbe  method  has  been  applied  on  an  extensive  domain,  encom¬ 
passing  123°-129°  E  and  26°-32°  N,  so  that  the  drifter  initial  loca¬ 
tions  are  centered  in  the  computational  domain  and  their 
evolution  during  the  experiment  is  captured.  Even  though  the  spe¬ 


cific  application  in  this  paper  is  limited  to  smaller  space  scales  and 
short  time  scales  of  1  or  2  days,  the  method  has  actually  been  ap¬ 
plied  to  a  longer  period,  extending  25  days  from  October  8  until 
November  2,  2007,  in  order  to  check  the  variability  of  the  process. 
A  total  of  27  of  the  30  deployed  drifters  have  been  used  in  the 
reconstruction,  since  the  other  three  do  not  last  through  the  whole 
analysis  period. 

The  input  data  for  the  reconstruction  are  the  two-dimensional 
velocity  fields  from  the  10-m  level  of  EAS-16  for  the  25  days  with 
At  =  1  h,  and  the  positions  of  the  27  drifters.  In  order  to  minimize 
the  errors  that  may  be  caused  by  the  time  difference  between  the 
two  input  data  sets,  the  drifter  positions  are  interpolated  to  the 
time  steps  of  the  velocity  fields.  The  three-dimensional  EAS-16 
data  set  is  then  used  to  statistically  project  the  corrected  velocity 
down  to  the  150  m  depth,  corresponding  to  the  ten  vertical  analy¬ 
sis  levels  of  the  model. 

The  basic  concept  of  LAVA  consists  of  correcting  the  model 
velocity  by  minimizing  the  difference  between  tbe  observed  drifter 
positions  and  the  model  forecasted  positions  from  numerical  tra¬ 
jectories  (see  Appendix  A).  The  local  correction  along  the  drifter 
path,  Au,  is  assumed  to  be  characterized  by  some  basic  correlation 
scales  in  time  (T)  and  space  (R),  characterizing  the  persistence  of 
the  flow.  Au  is  extended  with  scale  R  along  the  path  and  is  as¬ 
sumed  to  be  constant  in  time  over  T,  where  T  must  be  smaller  than 
or  of  the  same  order  as  the  Eulerian  time  scale  T^.  The  correction  is 
performed  over  sequences  with  steps  t  C  T,  that  are  also  required 
to  be  no  longer  than  the  Lagrangian  time  scale  T^  (Molcard  et  al., 
2003). 

These  main  scales  have  been  identified  by  computing  Eulerian 
and  Lagrangian  autocovariances  in  a  preliminary  analysis  of  the 
drifter  data  and  the  model  velocity  timeseries.  As  discussed  in 
Section  2,  the  region  of  interest  is  characterized  by  two  very  dif¬ 
ferent  regimes.  The  shelf  region  is  shallow  and  dominated  by 
weak  currents,  while  the  deep  area  is  dominated  by  the  presence 
of  the  Kuroshio  with  strong  and  highly  correlated  currents.  De¬ 
spite  these  differences,  the  time  scales  appear  relatively  homoge¬ 
neous  in  the  domain  with  dominant  mesoscale  Eulerian  time 
scales  Tf  on  the  order  of  1  week,  while  the  Lagrangian  time  scale 
Ti »  1-2  days,  which  is  typical  for  the  upper  ocean.  In  addition  to 
the  mesoscale,  there  is  also  a  strong  tidal  signal  with  a  period  of 
12  h.  In  the  present  application,  as  in  previous  work  (Taillandier 
et  al.,  2006a,  2008),  we  focus  on  corrections  at  the  mesoscale, 
since  higher  frequency  corrections  would  require  an  unrealisti¬ 
cally  high  data  coverage.  To  this  end,  the  correction  is  computed 
considering  the  12-h  low-pass  filtered  velocity  field.  A  number  of 
preliminary  tests  have  been  performed  considering  various 
parameter  values  for  T  and  t,  and  the  results  appear  to  be  quite 
robust.  The  results  reported  in  the  following  are  obtained  for 
T  =  6  h  and  3  h  for  the  deep  and  shallow  areas,  respectively.  T  is 
chosen  to  equal  t,  (T  =  x),  which  implies  that  the  procedure  is  ap¬ 
plied  over  sequences  of  1  step  only. 

As  for  the  space  scales,  we  expect  that  the  differences  between 
tbe  two  regimes  will  be  more  significant  and  that  they  are  sepa¬ 
rated  by  well  defined  frontal  dynamics.  Analysis  of  the  EAS-16 
velocity  field  suggests  that  the  300  m  isobath  approximates  the 
boundaiy  between  the  two  regimes  (Fig.  5(a)).  We  select  two  dif¬ 
ferent  space  scales,  20  km  for  the  deep  region  and  10  km  for  the 
shallow  one  (cf.  Taillandier  et  al.,  2008).  Notice  that  since  the  spa¬ 
tial  correlation  structure  is  modeled  as  the  solution  of  a  diffusion 
equation  (Berber  and  Rosati,  1989;  Weaver  and  Courtier,  2001), 
anisotropic  correlations  can  be  taken  into  account,  which  are  ex¬ 
pected  to  occur  for  strongly  sheared  currents  flowing  along 
topography. 

Fig.  5  gives  an  example  of  the  velocity  field  in  the  area  covered 
by  tbe  drifters.  Panel  (a)  shows  the  original  EAS-16  field,  panel  (b) 
tbe  corrected  field,  and  panel  (c)  the  difference  field.  The  correction 


184 


Y.  Chang  et  al./ Ocean  Modelling  36  (2011)  179-197 


0  IT  T  T  T  T  >r  ^  ^  ^  ^  ^  ^  ^  i  s>n^  ^  i  ^  .  - -^ 

126  126.2  126.4  126.6  126.8  127  127.2  127.4  127.6  127.8  128 


(b)  LAVA  velocity  field,  22:00  10/8/07  ~  21 :00  10/11/07 


(c)  Velocity  difference  (LAVA  -  EAS),  22:00  10/8/07  ~  21 :00  10/11/07 


Fig.  5.  Velocity  vectors  that  show  flow  patterns  in  the  region.  The  velocities  are  averaged  over  3  days  from  22:00  October  8  through  21:00  October  11, 2007.  (a)  Original  EAS- 
16  velocity;  (h)  corrected  LAVA  velocity;  (c)  velocity  difference  (LAVA  -  EAS-16).  Blue  lines  indicate  isobaths,  with  labels  in  m.  (For  interpretation  of  the  references  to  colour 
in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


induced  by  the  drifters  tends  to  decrease  the  velocity  in  the  shelf 
area,  while  it  increases  it  in  the  Kuroshio  and  in  the  transition  area. 
Once  we  have  computed  the  estimated  velocity  field  at  the  drif¬ 
ter  level,  we  need  to  calculate  its  vertical  projections  in  order  to 
estimate  the  velocities  at  other  depths.  To  do  this,  the  vertical  cor¬ 
relation  function,  C,  can  be  calculated  for  the  two  velocity  compo¬ 
nents  separately  as 


{u'(x,y,Zo,t)u’(x,y,z,t)) 

(u'(x,y,z„,t)u'(x,y,Zo,t))’ 


(1) 


where  u'  =  u  -  (u),  (u)  is  the  time  averaged  velocity  computed  over 
the  25  days  of  the  analysis  window,  and  Zg  is  the  depth  level  at 
which  the  velocity  has  been  estimated  by  the  reconstruction  pro¬ 
cess,  The  time  average  is  chosen  to  correspond  to  the  entire  analysis 
period  since  the  stratification  is  approximately  stationary  in  a  sta¬ 
tistical  sense. 

To  compute  C  in  Eq.  (1),  we  concentrate  on  the  area  where  the 
velocity  difference  between  the  EAS-16  and  the  LAVA  estimations 
is  significant,  i,e„  in  the  region  of  influence  of  the  drifters 
(Fig,  6(a)),  Since  the  water  depths  vary  inside  this  region,  we  divide 
the  region  into  six  sub-regions  according  to  the  water  depth 
(Fig,  6(b)).  Then  the  function  C  is  horizontally  averaged  inside  each 


sub-region  to  find  six  vertical  profiles  (Molcard  et  al.,  2005)  as 
shown  in  Fig.  6(c)  and  (d). 

The  vertical  correlation  is  then  used  to  estimate  the  vertical 
projections  of  velocity  corrections,  Au, 

Au(x,y,z,f)  =  Au(x,y,z„,t)  •  Q„(x,y,z),  (2) 

where  Au(x,y,Zo,t)  =  Ue(x,y,z„,t)  -  Umod(x.y.Zo.t).  and  Ug{x,y,Zo,t)  is 
the  velocity  estimated  by  LAVA,  while  Umoa(x,y,Zo,t)  is  the  EAS-16 
velocity. 

The  vertical  projection  of  the  estimated  velocity  field  is  then 
found  from 

Ue(x,y,z,t)  =Umod(x,y,z,f)  +  Au(x,y,z,t).  (3) 

In  Fig.  7,  one  example  of  the  vertically  projected  estimated  velocity 
corrections  at  150  m  is  compared  to  the  corresponding  velocity  cor¬ 
rections  at  the  drifter  level.  The  general  pattern  is  similar  between 
the  two  depths,  while  the  strength  is  weaker  at  the  deeper  level. 

3.2.  Constrained  normal  mode  analysis  (NMA) 

The  Normal  Mode  Analysis  (NMA)  data  blending  technique 
(Eremeev  et  al.,  1992a)  has  been  used  to  study  flow  in  Lake  Ontario 
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(a)  Velocity  difference  between  EAS  and  LAVA 


Longitude 


Fig.  6.  (a)  Velocity  vectors  that  show  the  difference  between  EAS-16  and  LAVA  averaged  over  25  days  over  bathymetry.  The  contours  show  the  water  depth  ranging  from 
500  m  (blue)  to  3000  m  (red),  (b)  Depth  contours  inside  the  region  where  the  velocity  difference  is  significant,  defining  six  distinct  sub-regions  according  to  water  depth  h. 
Region  1:  h  <  80  m.  Region  2:  80  c  h  <  100  m.  Region  3:  100  ^  h  <  150  m,  Region  4:  150  <  h  <  400  m,  Region  5:  400  ^h<  800  m.  Region  6:  h  >  800  m.  (c)  Vertical  correlation 
profiles  for  the  zonal  velocity  component  u.  (d)  Vertical  correlation  profiles  for  the  meridional  velocity  component  v.  Colors  in  the  correlation  plots  correspond  to  the  sub- 
regions  defined  in  (b).  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


(Rao  and  Schwab,  1981),  circulation  in  the  Black  Sea  (Eremeev 
et  al.,  1992b),  surface  circulation  in  Monterey  Bay  (Lipphardt 
et  al.,  2000),  among  many  other  applications.  The  method  has  also 
been  adapted  to  use  drifter  velocities  to  enhance  near-surface 
model  velocities  in  the  Gulf  of  Mexico  (Toner  et  al.,  2001).  Here, 
we  follow  Toner  et  al.  (2001),  using  velocities  derived  from  drifter 
trajectories  to  constrain  NMA  objective  mappings  of  model  veloc¬ 
ities  at  one  depth.  An  overview  of  the  methodology  and  the  param¬ 
eter  choices  is  provided  here,  while  details  are  given  in  Appendix  B. 

The  velocity  field  u  is  decomposed  into  the  sum  of  a  boundary 
solution  Utidry,  a  set  of  Dirichlet  modes  u°,  and  a  set  of  Neumann 
modes  uJJ,  (see  details  in  Appendix  B).  Dirichlet  modes  are  eigen¬ 
functions  of  the  Helmholtz  equation  with  Dirichlet  boundary  con¬ 
ditions,  while  Neumann  modes  are  eigenfunctions  of  the 
Helmholtz  equation  with  Neumann  boundary  conditions.  Thus, 
the  velocity  estimate  is  given  by 

No 

Ue{x,y,  t)  =  UMry(x,y,  t)  +  •  •  •  a„(tX(x,y) 

n=l 

+  ^h„(tX(x,y),  (4) 

m=l 

where  No  is  the  number  of  Dirichlet  modes  and  Njv  is  the  number  of 
Neumann  modes,  and  and  are  the  coefficients.  The  boundary 
solution  is  determined  separately  by  matching  normal  flow  on  the 
boundary  and  solving  for  the  velocity  potential. 


The  coefficients  are  determined  as  the  solution  to  a  least 
squares  problem  minimizing 

||tle  ~  ttmodll ;  (5) 

where  ||  ■  |[  denotes  the  2-norm.  (This  is  equivalent  to  the  minimiza¬ 
tion  problem  (B.IO).) 

NMA  has  several  useful  properties;  it  reduces  the  number  of  de¬ 
grees  of  freedom,  identifies  the  dominant  length  scales  in  a  dy¬ 
namic  region,  and  applies  systematic  spatial  smoothing  through 
the  choice  of  a  length  scale  cut-off  for  the  eigenfunctions.  It  can 
also  be  applied  to  irregularly  spaced  data  and  data  from  disparate 
sources,  which  is  of  primary  importance  in  our  context.  There  are 
various  choices  for  the  introduction  of  the  obseiA/ed  velocities.  At 
one  end  of  the  spectrum,  they  can  be  treated  identically  to  the 
model  velocities,  introducing  additional  components  to  the  cost 
function.  In  many  cases,  however,  observations  are  taken  to  be 
more  certain,  and  this  can  be  accounted  for  in  the  set-up  of  the 
minimization  problem  through  differential  weighting.  At  the  other 
end  of  the  spectrum— the  method  employed  here— the  obseiT/a- 
tions  are  imposed  as  hard  constraints  on  the  cost  function. 

The  primary  parameter  that  enters  the  set-up  of  this  type  of 
NMA  is  the  choice  of  the  number  of  modes  or  the  minimum  re¬ 
solved  length  scale.  If  too  few  modes  are  used,  the  solution  may 
be  forced  to  be  smoother  than  is  realistic.  At  the  same  time,  the 
hard  constraints  may  result  in  undesirable  large-scale  spatial  oscil¬ 
lations.  (In  the  extreme  case,  the  constrained  minimization  has  no 
solution.)  On  the  other  hand,  if  too  many  modes  are  used,  then,  in 
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(a)  velocity  difference  at  z=10  m,  Oct.  12,  06:00 


Fig.  7.  Velocity  difference  between  EAS-16  and  LAVA  at  (a)  the  drifter  level  and  (b)  the  vertical  projection  at  z=  150  m,  over  bathymetry.  The  contours  show  the  water  depth 
ranging  from  500  m  (blue)  to  3000  m  (red).  This  example  is  taken  at  6:00  Oct.  12,  2007.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is 
referred  to  the  web  version  of  this  article.) 


effect,  the  model  is  weighted  more  heavily  and  may  reduce  the 
sphere  of  influence  of  each  obsen/ation  too  much.  Choosing  the 
correct  length  scale  cut-off  to  minimize  spurious  long  distance  ef¬ 
fects  yet  maximize  meaningful  impact  of  the  observations  is  some¬ 
thing  of  an  inexact  art,  usually  involving  some  trial  and  error. 

Here  we  have  chosen  the  length  scale  cut-off  to  be  at  0.125°,  the 
equivalent  of  twice  the  grid  spacing,  and  the  domain  125.5°-128.5° 
E  by  27.75°-30.75°  N.  This  results  in  No  =  576  Dirichlet  and 
Nn  =  624  Neumann  modes.  Given  such  a  rich  data  set  of  observa¬ 
tions  as  the  LWAD07  SVP  drifters,  we  actually  have  on  occasion 
multiple  observations  that  are  closer  than  the  smallest  resolved 
spatial  scale.  While  the  machineiy  of  NMA  is  capable  of  solving 
the  constrained  least  squares  minimization  nonetheless— due  to 
the  large  number  of  modes  relative  to  the  number  of  observa¬ 
tions— the  result  is  not  necessarily  sensible.  Consequently,  we 
implement  a  filter  on  the  observations,  imposing  0.125°  as  the 
minimum  distance  between  included  data  points. 

Theoretically,  it  is  possible  to  extend  NMA  into  three  dimen¬ 
sions  and  carry  out  a  three-dimensional  constrained  minimization. 
However,  this  is  computationally  costly.  As  all  observations  in  our 
case  are  restricted  to  a  single  depth,  we  apply  a  constrained  NMA 
to  the  nearest  model  level  (10  m)  and  project  the  innovations 
down  (and  up)  in  the  water  column  using  linear  regression.  At  each 
time,  an  unconstrained  NMA  map  is  found  for  all  levels,  in  addition 
to  the  constrained  NMA  map  for  the  reference  level.  For  each  non¬ 


reference  level  separately,  a  linear  regression  coefficient  C  is  deter¬ 
mined  for  each  of  the  amplitudes  for  each  24  h  window,  relative  to 
the  corresponding  reference  level  amplitudes  (Eq.  (6)).  This  coeffi¬ 
cient  is  used  to  project  the  innovation  increment  from  the  refer¬ 
ence  level  to  the  other  level,  which  is  then  added  to  the 
unconstrained  NMA  map  of  the  non-reference  level  (Eq.  (7)).  Math¬ 
ematically,  this  process  can  be  summarized  as  follows,  where  a 
represents  an  amplitude  time  series,  k  indexes  the  amplitudes,  Zg 
refers  to  the  reference  level,  subscript  e  marks  amplitude  estimates 
from  the  constrained  NMA,  primes  indicate  perturbations  from  a 
24-h  mean,  and  superscript  T  denotes  the  matrix  transpose: 


Cz.(/GZ) 
ae(/c,z) : 


a'(k,zya'(k,Zo) 
a'{k,ZoYa'{k,Zo)' 
a{k,z)  +  Q„(k,z)(ae{k, 


Zo) 


a(k,Zo)). 


(6) 

(7) 


This  method  of  statistical  projection  is  identical  to  that  described  in 
Section  3.1,  with  two  exceptions:  (1)  The  correlation  is  carried  out 
on  NMA  amplitudes  instead  of  directly  on  the  velocities,  and  (2) 
the  time  average  is  taken  over  24  h  instead  of  over  25  days. 

Note  that  the  boundary  solution  of  neither  the  reference  level 
nor  any  other  level  Is  affected  by  the  observations.  This  is  a  limita¬ 
tion  of  the  NMA  method,  but  does  not  cause  too  much  trouble  as 
long  as  the  boundary  of  the  domain  used  is  far  enough  away  from 
the  area  of  interest,  a  typical  restriction  to  avoid  unintended 
boundary  effects  (for  more  details,  please  refer  to  Appendix  B). 
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The  method  has  been  applied  to  create  hourly  NMA  velocity 
fields  for  the  time  period  from  12:00  October  8,  2007  to  11:00 
October  15,  2007.  The  first  observation  is  available  for  22:00 
October  8,  2007,  from  which  point  onwards  the  constrained  NMA 
algorithm  is  used.  The  number  of  obsereations  blended  each  hour 
varies  from  1  to  14.  Recall  that  spatial  filtering  of  observations  was 
applied  to  avoid  imposing  unresolved  spatial  variability.  Conse¬ 
quently,  the  full  set  of  thirty  drifters  was  never  used,  even  during 
the  period  when  all  were  reporting. 

The  chosen  least  resolved  spatial  scale  is  fairly  small,  so  that  the 
impact  of  the  observations  remains  somewhat  localized.  This  is,  of 
course,  particularly  true  when  few  observations  are  absorbed. 
Fig.  8  shows  the  velocity  fields  from  the  original  EAS-16  model  and 
from  the  constrained  NMA,  as  well  as  the  differences  in  the  fields, 
for  two  sample  days  with  5  and  13  observations  used,  respectively. 
Fig.  9  illustrates  the  effects  on  a  deeper  level,  in  this  case  the  one  at 
40  m.  We  can  see  that  at  the  lower  levels,  too,  the  corrections  are  re¬ 
stricted  to  the  geographic  area  immediately  near  the  observations 
but  are  smaller  than  at  the  reference  level.  Both  of  these  characteris¬ 
tics  are  desirable,  since  velocities  at  depth  tend  to  be  smaller. 

3.3.  Internal  consistency  tests  using  drifter  trajectories 

A  first  test  on  the  corrected  velocity  fields  obtained  with  the  two 
methods  in  Sections  3.1  and  3.2  is  performed  comparing  synthetic 


trajectories  with  the  observed  in  situ  drifter  trajectories.  Since  the 
velocity  correction  is  based  on  the  drifter  data  themselves,  the  test 
is  mostly  aimed  at  verifying  internal  consistency,  and  a  significant 
improvement  is  expected  to  occur  for  the  trajectories  in  the  cor¬ 
rected  fields  with  respect  to  the  original  EAS-16  ones. 

The  synthetic  trajectories  are  numerically  advected  using  the 
model  velocity  fields  starting  from  the  launch  locations  of  the 
in  situ  drifters.  Fig.  10  shows  two-day  trajectories  after  the 
launches  based  on  the  following  fields:  EAS-16  (red,  left  panel), 
NMA  corrected  (blue,  center  panel),  and  LAVA  corrected  (green, 
right  panel).  In  all  the  panels,  the  synthetic  trajectories  are  shown 
together  with  the  in  situ  ones  (black).  It  is  evident  from  visual 
inspection  that  both  methods  lead  to  a  substantial  improvement. 
The  NMA  and  LAVA  trajectories  appear  significantly  closer  to  the 
observations  than  the  EAS-16  ones,  capturing  the  two  different 
flow  regimes,  on  the  shelf  and  in  the  Kuroshio,  respectively.  The 
corrected  trajectories  on  the  shelf  are  characterized  by  tidal  motion 
and  reduced  drift,  while  those  caught  in  the  Kuroshio  exhibit  a 
swift  northeastern  drift,  similar  to  the  observations. 

This  result  is  quantified  by  computing  the  rms  error  in  the  mod¬ 
eled  versus  observed  positions,  averaged  over  all  the  available  drift¬ 
ers,  as  a  function  of  time: 


EAS-16 


Const! .  NMA 


Difference 


Fig.  8.  Comparison  of  velocity  fields  at  the  10  m  reference  level  from  (left)  EAS-16  and  (middle)  constrained  NMA.  The  panels  on  the  right  depict  the  differences  (NMA  -  EAS- 
1 6);  note  the  different  vector  scale  used  there.  Red  arrows  show  the  assimilated  observations.  The  top  three  panels  are  for  1 2:00  October  9, 2007,  when  only  five  observations 
were  used,  while  the  lower  three  panels  are  for  12:00  October  12,  2007,  when  13  obsei'vations  were  used.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  article.) 


188 


Y.  Chang  et  al./ Ocean  Modelling  36  (2011)  179-197 


EAS-16  Proj.  Constr.  NMA  Difference 


Fig.  9.  Comparison  of  velocity  fields  at  the  40  m  level  from  (left)  EAS-16  and  (middle)  projected  constrained  NMA.  The  panels  on  the  right  depict  the  differences  (NMA  -  EAS- 
16);  note  the  different  vector  scale  used  there.  Red  arrows  show  the  assimilated  observations.  These  plots  are  for  12:00  October  12. 2007,  like  the  bottom  panels  in  Fig.  8.  (For 
interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


Fig.  10.  Modeled  drifter  trajectories  using  the  EAS-16  (red.  a).  NMA  corrected  (blue,  b),  and  LAVA  corrected  (green,  c)  model  velocity  fields,  compared  to  in  situ  data  (black) 
over  two  days  following  launches.  Purple  circles  mark  initial  locations.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 


where  robs(n)  and  TmodW  are  the  positions  of  the  nth  observed  and 
modeled  trajectories,  respectiveiy,  at  the  set  time  and  Ndd/t  is  the 
number  of  drifters  availabie  at  that  time.  The  error  for  the  three  sets 
of  trajectories  computed  from  EAS-16,  NMA,  and  LAVA,  respec¬ 
tively,  are  shown  in  Fig.  11,  where  the  standard  deviation  a  is  aiso 
piotted  to  indicate  the  range  of  the  trajectory  errors.  While  the  EAS- 
1 6  error  (red  line)  increases  aimost  iineariy,  reaching  approximateiy 
21  km  after  24  h,  the  NMA  and  LAVA  curves  (blue  and  green, 
respectively)  show  a  slower  rate  of  growth,  reaching  approximateiy 
9.5  km  for  NMA  and  5.5  km  for  LAVA  after  24  h. 

Fig.  12  displays  the  separation  distance  as  a  function  of  time  for 
the  individual  drifters  between  the  observed  trajectory  and  that 
modeled  by  each  of  the  methods.Note  that  the  distributions  across 
drifters  are  skewed  with  a  tail  at  the  upper  end.  Separation  dis¬ 
tances  for  each  of  the  improved  models  are  more  clustered  towards 
the  lower  end.While  NMA  has  more  trajectories  that  stay  within 
3  km  of  the  observations  after  24  h  (18)  than  LAVA  (10),  it  also  re¬ 
sults  in  larger  excursions  with  a  maximum  of  27.5  km,  versus 
17.1  km  for  LAVA. 

We  characterize  the  improvement  or  gain  due  to  the  corrected 
veiocities  with  a  simpie  bulk  metric  C,  given  by  the  normalized  dif¬ 
ference  between  the  rms  errors  at  24  h  for  the  EAS-16  and  the  cor¬ 
rected  models.  The  24  h  horizon  is  chosen  as  particuiarly  reievant 


time,  h 

Fig.  11.  Comparison  of  rms  errors  for  drifter  trajectories  modeled  using  EAS-16 
(red),  NMA  corrected  (blue)  and  LAVA  corrected  (green)  velocity  fields.  Time  is 
measured  from  each  drifter's  launch.  Error  bars  show  ±  one  standard  deviation.  (Eor 
interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  article.) 

for  operational  applications.  This  provides  the  foilowing  values  for 
NMA  and  LAVA,  respectively: 


y.  Chang  et  al.  /  Ocean  Modelling  36  (2011)  179-197 


189 


Gnma  —  1  ~  Enma/Eeas  ~  55%,  (9) 

Clava  ~  ^  —  Elava/Eeas  ~  75%.  (10) 

4.  Assessment  of  velocity  reconstruction  using  sonobuoy 
trajectories 

The  corrected  velocity  fields  obtained  using  SVP  drifter  data  as 
described  in  Section  3  are  quantitatively  tested  here  against  data 
from  the  CDMR  sonobuoy  trajectories.  The  testing  methodology 
is  similar  to  the  one  used  in  Section  3.3,  i.e.  synthetic  trajectories 
are  computed  based  on  the  velocity  fields  and  compared  with 
the  observed  trajectories.  There  is,  however,  a  conceptual  differ¬ 
ence  between  the  present  analysis  and  the  one  performed  in  Sec¬ 
tion  3.3  using  drifters,  since  the  sonobuoy  trajectories  provide  an 
independent  data  set,  that  has  not  been  used  in  the  velocity  recon¬ 
struction.  We  recall  that  sonobuoy  trajectories  differ  significantly 
from  the  drifter  ones  (see  Fig.  3)  due  to  different  depth  dependent 
drag  and  different  initial  conditions.  The  comparison  performed 
here  therefore  provides  an  independent  quantitative  test  of  the 
quality  of  the  reconstructed  velocities. 

Synthetic  sonobuoy  trajectories  are  calculated  using  the  Sono¬ 
buoy  Field  Drift  Model  (SFDM),  see  Hammond  (2005,  2008a,b) 
and  Appendix  C  for  details.  The  code  was  developed  as  a  Navy  mis¬ 
sion  planning  tool  under  the  ODDAS  (Optimal  Deployment  of  Drift¬ 
ing  Acoustic  Sensors)  program.  SFDM  models  the  hydrodynamic 
response  of  sonobuoys— represented  as  a  distributed  series  of  sur¬ 
face  and  subsurface  drag/mass  bodies  and  cable— to  a  spatially  and 
temporally  varying  water  velocity  field.  While  aerodynamic  drag 
on  the  surface  buoy  is  an  integral  part  of  the  model,  no  wind  drag 
was  taken  into  account  here,  since  wind  field  data  was  not  avail¬ 


able.  However,  over  90%  of  the  drag  area  of  a  typical  sonobuoy  lies 
below  the  surface  of  the  water,  and  previous  studies  (Goughian, 
1975)  have  shown  that  wind  typically  contributes  less  than  2%  to 
the  sonobuoy  drift. 

SFDM  simulations  were  conducted  for  the  LWAD07  sonobuoy 
deployments  using  the  three  model  velocity  fields;  EAS-16,  NMA 
corrected,  and  LAVA  corrected.  These  modeled  trajectories  in  com¬ 
parison  to  those  from  the  observed  sonobuoy  data  are  plotted  in 
Fig.  13.  Qualitatively  the  SFDM  modeled  results  using  the  corrected 
velocity  fields  show  better  correlation  to  the  in  situ  data.  As  for  the 
drifter  data  in  Fig.  10,  the  EAS-16  predictions  tend  to  be  fairly  uni¬ 
form  in  space,  whereas  both  of  the  corrected  velocity  fields  are  bet¬ 
ter  able  to  predict  the  looping  behavior  of  the  sonobuoys  deployed 
in  the  northwestern  corner,  while  also  increasing  the  energy  of  the 
trajectories  in  the  southwestern  corner  and  in  the  Kuroshio  (even 
though  the  drift  of  the  current  is  still  underpredicted). 

The  separation  distance  A  between  the  in  situ  data  and  the 
SEDM  results  is  computed  for  all  trajectories  and  plotted  as  a 
function  of  time  in  Eig.  14  for  the  entire  operating  life  of  each  sono¬ 
buoy.  This  data  supports  the  qualitative  conclusion  that  the  cor¬ 
rected  velocity  fields  improve  the  sonobuoy  drift  predictions. 
While  the  majority  of  the  sonobuoy  trajectories  predicted  using 
the  corrected  velocity  fields  (blue  and  green)  exhibit  separation 
distances  less  than  20  km  over  the  whole  life  span,  there  are  veiy 
few  trajectories  with  distances  consistently  less  than  20  km  when 
the  EAS-16  output  (red)  Is  used  for  SEDM  predictions. 

In  Elg.  15,  the  rms  separation  distance  averaged  over  all  trajec¬ 
tories,  corresponding  to  the  error  E  defined  in  Eq.  (8),  is  shown  for 
each  of  the  model  velocity  fields.  Notice  that  £(f)  tends  to  plateau 
for  all  velocities  after  approximately  20  h,  at  which  point  a  few 
trajectories  with  high  errors  die,  as  shown  by  the  individual  A(t) 


EAS-16 


0  3  6  9  12  15  18  21  24 

time  from  launch  (h) 

Fig.  12.  Timeseries  of  separation  distances  between  in  situ  and  modeled  drifter  trajectories  using  EAS-16  (red),  NMA  corrected  (blue),  and  LAVA  corrected  (green)  velocity 
fields.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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(a)  EAS  vs.  sonobuoy  GPS 


(b)  NMA  vs.  sonobuoy  GPS 


(c)  LAVA  vs.  sonobuoy  GPS 
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Fig.  13.  Modeled  (SFDM)  sonobuoy  trajectories  using  the  EAS-16  (red,  a).  NMA  corrected  (blue,  b)  and  LAVA  corrected  (green,  c)  model  velocity  fields,  compared  to  in  situ 
data  (black)  for  two  days  after  deployment.  Purple  circles  mark  initial  locations.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to 
the  web  version  of  this  article.) 


in  Fig.  14.  Both  the  NMA  and  LAVA  corrected  velocity  fields  show  a 
significantly  improved  performance  over  the  EAS-16  results. 

The  gain  at  24  h  is  computed  from  E  as  in  Section  3.3,  giving  for 
NMA  and  LAVA,  respectively, 

25%,  (11) 

Gmva  ~  50%.  (12) 

These  values  are  smaller  than  those  obtained  in  Section  3.3  for  drift¬ 
ers,  as  is  expected,  given  that  the  sonobuoy  data  are  not  used  in  the 
reconstruction.  The  decrease  is  relatively  small,  though,  on  the  or¬ 
der  of  25%  for  both  fields,  suggesting  that  the  results  of  the  methods 
are  robust  and  provide  useful  results  not  only  at  the  drifter  level  but 
also  throughout  the  water  column. 

Another  practical  metric  to  characterize  the  performance  is  the 
time  to  separation  Tj,  for  some  distance,  L,  between  the  modeled 


and  in  situ  trajectories.  For  the  purposes  of  this  study  a  distance 
of  roughly  2  grid  spacings  or  15  km  has  been  chosen  for  L  The 
T]5  metric  cannot  be  used  to  make  average  comparisons  of  overall 
field  performance,  since  in  many  cases  the  sonobuoy  reaches  the 
end  of  its  life  (EOL)  before  separating  more  than  15  km  from  the 
modeled  trajectories.  However,  it  can  be  used  to  compare  the  per¬ 
formance  for  individual  sonobuoys. 

It  is  instructive  to  examine  some  specific  examples  to  illustrate 
the  effectiveness  of  using  the  corrected  velocity  fields  for  sonobuoy 
drift  prediction.  The  following  three  examples  demonstrate  typical 
behavior  of  sonobuoys  deployed  in  three  regions  of  the  test  area. 
The  in  situ  data  and  SFDM  results  using  the  three  model  velocity 
fields  are  plotted  together  and  the  separation  distance  at  24  h, 
A24,  as  well  as  the  T15  metric  are  compared. 

Buoy  #1  is  shown  in  Fig.  16.  This  sonobuoy  was  deployed  in  the 
transition  region  between  the  two  regimes  of  the  shelf  and  of  the 


Fig.  14.  Timeseries  of  separation  distances  between  in  situ  and  modeled  (SFDM)  sonobuoy  trajectories  using  EAS-16  (red),  NMA  corrected  (blue),  and  LAVA  corrected  (green) 
velocity  fields.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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Fig.  15.  Comparison  of  rms  errors  for  modeled  sonobuoy  trajectories  (SFDM)  using 
EAS-16  (red),  NMA  corrected  (blue),  and  LAVA  corrected  (green)  velocity  fields. 
Time  is  measured  from  each  sonobuoy’s  launch.  Error  bars  show  ±  one  standard 
deviation.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the 
reader  is  referred  to  the  web  version  of  this  article.) 


Kuroshio.  Here  the  in  situ  data  show  that  the  sonobuoy  initially 
moved  to  the  southeast  before  getting  caught  up  in  the  Kuroshio 
and  drifting  rapidly  to  the  northeast.  The  drift  predictions  using 
the  EAS-16  velocity  field  data  do  not  capture  this  type  of  motion. 
The  24-h  separation  distance  is  almost  14  km  increasing  to 
41  km  by  the  end  of  the  buoy’s  life  (61.5  h).  However,  drift  predic¬ 
tions  made  using  either  corrected  velocity  field  do  exhibit  the  cor¬ 
rect  drift  pattern.  At  24  h  both  predictions  show  low  separation 
distances.  Drift  velocity  using  the  corrected  velocity  fields  is  gener¬ 
ally  lower  than  observed  once  the  sonobuoy  is  in  the  Kuroshio,  and 
separation  distance  increases  over  the  life  of  the  buoy.  The  LAVA 
predictions  show  better  performance  in  this  case  than  the  NMA 
predictions,  staying  within  15  km  of  the  observed  position  for  over 
51  h,  versus  about  41  h  using  NMA  corrected  velocity  fields  and 
26  h  for  EAS-16, 

Buoy  #8  (see  Fig,  17)  was  deployed  in  the  southeast  corner  of 
the  test  site  where  the  influence  of  the  Kuroshio  was  dominant. 
After  initially  drifting  to  the  ENE  at  a  high  velocity  for  roughly 
the  first  8  h,  Buoy  #8  slows  slightly  and  drifts  NNE  until  it  dies. 
None  of  the  SFDM  results  predict  the  initial  ENE  motion.  All  three 
data  sets  do  predict  that  Buoy  #8  will  drift  to  the  NNE,  although  at 
a  lower  velocity  than  the  obseiA/ed  data.  This  initial  divergence  is 
why  none  of  the  predicted  trajectories  were  able  to  remain  within 
15  km  of  the  observed  data  for  more  than  10,2  h.  Nonetheless,  due 
to  the  better  drift  speed  correlation,  the  corrected  velocity  predic¬ 
tions  show  significant  improvement  in  separation  distance  after 
24  h;  37.3  km  with  the  EAS-16  output  and  roughly  20  km  for  both 
corrected  velocity  fields. 

Buoy  #18  is  plotted  in  Fig.  18.  This  buoy  was  deployed  in  the 
northwest  corner  of  the  test  site.  As  was  typical  for  sonobuoys  de¬ 
ployed  in  that  region,  it  was  observed  to  loop  with  a  semi-diurnal 
period  and  slow  mean  translation  to  the  southwest.  The  SFDM  re¬ 
sults  using  the  EAS-16  velocity  field  fail  to  replicate  this  motion. 
Although  there  is  a  semi-diurnal  looping  present,  there  is  a  signif¬ 
icant  mean  translation  to  the  northeast  such  that  by  the  end  of  the 
buoy  life  (68  h)  the  EAS-16  predicted  position  is  over  30  km  from 
the  observed.  The  corrected  velocity  field  predictions  perform 
much  better,  as  they  capture  the  general  character  of  the  obseiT/ed 
motion.  Separation  distances  are  more  a  function  of  phase  differ¬ 
ences  in  the  elliptical  motion  between  the  predicted  and  obseiT/ed 
positions.  The  NMA  corrected  results  indicate  separations  dis¬ 
tances  less  than  1,5  km  for  the  life  of  the  buoy,  while  the  LAVA  cor¬ 
rected  results  remain  less  than  5  km  from  observations. 
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Fig.  16.  Plot  of  observed  sonobuoy  trajectory  (black)  and  SFDM  modeled  trajecto¬ 
ries  using  EAS-16  (red),  NMA  corrected  (blue),  and  LAVA  corrected  (green)  velocity 
fields  for  Buoy  #1.  Bathymetry  is  shown  in  shades  of  blue  (in  m).  (for  interpretation 
of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 
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Fig.  17.  Plot  of  observed  sonobuoy  trajectory  (black)  and  SFDM  modeled  trajecto¬ 
ries  using  EAS-16  (red),  NMA  corrected  (blue),  and  LAVA  corrected  (green)  velocity 
fields  for  Buoy  #8.  Bathymetry  is  shown  in  shades  of  blue  (in  m).  (for  interpretation 
of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 


The  performance  of  the  three  models  is  summarized  in  Table  1 
for  each  of  the  sonobuoys  deployed  during  LWAD07  in  terms  of 
separation  distance  after  24  h  and  at  the  end  of  the  sonobuoy  life. 


5.  Summary  and  concluding  remarks 

In  this  paper  we  provide  a  quantitative  assessment  of  the 
improvement  in  trajectory  prediction  obtained  with  velocity  fields 
corrected  using  Lagrangian  data.  We  consider  the  specific  application 
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Fig.  18.  Plot  of  observed  sonobuoy  trajectory  (black)  and  SFDM  modeled  trajecto¬ 
ries  using  EAS-16  (red),  NMA  corrected  (blue),  and  LAVA  corrected  (green)  velocity 
fields  for  Buoy  #18.  Bathymetry  is  shown  in  shades  of  blue  (in  m).  (For 
interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is 
referred  to  the  web  version  of  this  article.) 


of  predicting  sonobuoy  trajectories  launched  during  the  LWAD07 
experiment  off  Taiwan. 

The  velocity  fields  of  the  data  assimilating  model  EAS-16  are 
corrected  using  data  from  SVP  drifters  launched  during  the  same 
time  period  and  drogued  at  15  m.  The  corrections  are  statistically 
propagated  from  the  drifter  depth  throughout  the  water  column 
using  the  statistics  of  the  (uncorrected)  EAS-16  fields.  Synthetic 
sonobuoy  trajectories  are  then  computed  from  the  EAS-16  and 
from  the  corrected  velocity  fields  using  an  appropriate  drift  model 
that  describes  the  hydrodynamic  response  of  sonobuoys  in  terms 
of  distributed  drag  in  the  upper  27  m,  and  the  synthetic  trajectories 
are  quantitatively  compared  with  the  observed  ones.  The  region  of 
interest  is  especially  challenging  since  it  is  situated  at  the  edge  of 
the  Kuroshio.  Both  drifter  and  sonobuoy  data  sets  show  the  pres¬ 
ence  of  two  distinct  flow  regimes  occurring  in  a  restricted  domain 
(order  of  1/2°  square)  and  separated  by  a  sharp  transition  area.  The 
western  side  is  characterized  by  shelf  dynamics  with  a  significant 
tidal  component,  while  the  more  eastern  part  is  dominated  by 
the  highly  energetic  meandering  of  the  Kuroshio.  The  drifter  and 
sonobuoy  data  sets,  even  though  qualitatively  similar,  have  signif¬ 
icant  quantitative  differences,  with  the  sonobuoy  trajectories  being 
overall  slower  and  exhibiting  slightly  different  orientations.  This  is 
likely  due  to  the  fact  that  sonobuoys  respond  differently  to  cur¬ 
rents  at  different  depths  and  that  they  have  been  launched  at 
slightly  different  points  and  times  than  the  drifters.  The  EAS-16 
model  does  not  capture  the  presence  of  the  two  regimes,  and  its 
synthetic  trajectories  show  a  more  homogeneous  northeastward 
motion  with  a  strong  tidal  component.  This  is  not  surprising,  given 
the  high  nonlinearity  of  the  flow,  that  makes  its  detailed  structure 
difficult  to  predict,  and  the  small  size  of  the  target  area  with  re¬ 
spect  to  the  resolution  of  the  assimilation  observing  system. 

Two  different  methods  for  velocity  reconstruction  have  been 
considered.  One  method,  LAVA,  is  based  on  a  variational  approach. 
It  corrects  the  velocity  held  by  requiring  minimization  of  the  dis¬ 
tance  between  obseiA/ed  and  modeled  trajectories.  The  other  meth¬ 
od,  NMA,  is  based  on  an  expansion  into  normal  modes.  It  uses  the 
velocities  computed  along  the  trajectories  as  hard  constraints  on 
the  cost  function. 


The  results  are  qualitatively  compared  with  the  observations  by 
visual  inspection  of  the  trajectories.  A  quantitative  analysis  is  then 
performed,  considering  as  main  metric  the  rms  separation  (error) 
between  the  observed  and  synthetic  trajectories  after  24  h.  An¬ 
other  metric  is  also  considered,  namely  the  time  to  separation  be¬ 
tween  observed  and  synthetic  trajectories  of  15  km  (roughly 
corresponding  to  twice  the  model  horizontal  grid  spacing). 

We  first  test  for  internal  consistency  of  the  two  methods,  con¬ 
sidering  the  SVP  drifter  trajectories.  While  the  EAS-16  trajectories 
do  not  reproduce  the  two  dynamical  regimes,  the  trajectories 
based  on  the  corrected  fields  appear  quite  similar  to  the  observed 
ones  and  show  the  significant  difference  between  the  two  regimes. 
The  error  is  decreased,  with  a  gain  of  approximately  75%  for  LAVA 
and  55%  for  NMA.  An  improvement  is  indeed  to  be  expected,  since 
the  fields  are  corrected  based  on  the  drifters.  Nevertheless,  our 
analysis  provides  a  quantitative  test  that  the  methods  produce 
consistent  results. 

We  then  consider  the  independent  sonobuoy  data  set  and  per¬ 
form  a  similar  trajectoiy  comparison.  As  for  the  drifters,  both  LAVA 
and  NMA  based  trajectories  are  qualitatively  much  more  similar  to 
the  observed  ones  than  the  EAS-16  based  ones.  Quantitatively,  the 
error  decreases  less  than  for  the  drifters,  but  still  significantly: 
approximately  50%  for  LAVA  and  25%  for  NMA.  For  the  separation 
time  metric,  a  bulk  value  for  the  whole  data  set  cannot  be  com¬ 
puted  since  many  trajectories  do  not  reach  a  separation  of  15  km 
during  their  entire  life  (which  is  on  the  order  of  two  days).  Compar¬ 
isons  are  performed  considering  single  trajectories,  and  a  signifi¬ 
cant  general  improvement  is  found  for  both  LAVA  and  NMA  with 
respect  to  EAS-16.  Best  results  appear  to  be  obtained  for  both 
methods  for  trajectories  situated  on  the  shelf,  whose  velocities 
tend  to  be  overestimated  by  the  model.  The  drifter-based  correc¬ 
tion  reduces  the  average  velocity  on  the  shelf,  and  the  resulting 
trajectories  tend  to  loop  almost  in  place,  in  keeping  with  the  obser¬ 
vations.  Trajectories  in  the  Kuroshio  are  on  average  more  energetic 
thanks  to  the  correction,  although  in  some  cases  not  energetic 
enough  compared  to  the  observed  ones.  The  correction  in  the 
Kuroshio  region  might  be  less  effective  due  to  the  fact  that  drifters 
there  move  very  quickly,  resulting  in  sparser  coverage,  not  always 
adequate  for  the  correction. 

Both  methods  produced  significant  improvement  in  recon¬ 
structing  the  sonobuoy  trajectories.  Even  though  based  on  two 
completely  different  approaches,  they  provide  qualitatively  similar 
results.  Quantitatively,  LAVA  provides  a  greater  error  reduction. 
This  is  likely  due  to  at  least  two  factors.  Results  from  previous  work 
(Molcard  et  al.,  2003)  show  that  trajectory-based  corrections  are 
generally  more  effective  than  corrections  relying  on  derived  veloc¬ 
ities,  at  least  when  the  time  steps  of  correction  are  a  sizable  frac¬ 
tion  of  the  Lagrangian  time  scales  Q,  as  in  this  case.  Conceptually 
this  is  due  to  the  fact  that  Eulerian  and  Lagrangian  velocities  are 
different  at  those  scales  and  the  Lagrangian  observational  operator 
is  more  appropriate,  allowing  the  minimization  of  the  difference 
between  observed  and  modeled  trajectories. 

A  more  specific  difference  between  the  methods  is  that  the 
NMA  methodology  used  here  imposes  a  hard  constraint,  requiring 
that  the  reconstruction  exactly  matches  the  drifter  obsei'vations, 
while  LAVA  optimally  weighs  the  observations  with  no  such 
requirement.  Hence  LAVA  allows  for  some  error  in  the  observa¬ 
tions  whereas  NMA  does  not.  A  further  consequence  of  this  partic¬ 
ular  formulation  of  NMA  is  that  not  all  the  available  data  is 
ingested:  At  most  14  of  the  30  drifters  were  used  at  any  one 
time.  In  addition,  NMA  results  in  a  more  localized  correction  than 
LAVA.  Of  course,  NMA  can  be  readily  adapted  to  permit  observa¬ 
tional  error  by  converting  the  constrained  optimization  into  a 
weighted  optimization  problem.  This  would  likely  improve  its 
performance.  LAVA,  moreover,  could  probably  be  further  opti¬ 
mized  as  well.  The  present  purpose,  however,  was  to  demonstrate 
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Table  1 

Summary  of  separation  distance  between  modeled  and  observed  sonobuoy  trajectories  after  24  h  and  at  end  of  buoy  life  for  all  three  velocity  fields.  (Note  that  Buoy  #  22  never 
deployed  and  is  hence  not  listed  below.) 


Buoy  Lifespan  (h)  Separation  distance  after  24  h  (km)  Separation  distance  at  EOL  (km) 


EAS-16 

NMA 

LAVA 

EAS-16 

NMA 

LAVA 

1 1 

61.50 

13.8 

4.5 

3.6 

41.4 

37.5 

21.3 

2 

24.50 

6.1 

5.8 

2.5 

6.3 

6.0 

2.4 

3 

22.00 

n/a 

n/a 

n/a 

36.3 

28.3 

9.4 

4 

27.25 

11.5 

5.6 

4.4 

14.2 

7.3 

4.  6 

5 

26.00 

13.5 

7.6 

2.  9 

15.2 

8.6 

3.0 

6 

35.75 

13.9 

16.4 

18.7 

25.9 

24.7 

19.8 

7 

41.25 

9.4 

0.4 

2.5 

14.8 

1.8 

3.  8 

1  8 

38.75 

37.3 

20.8 

20.1 

40.9 

21.7 

20.8 

9 

68.00 

13.9 

14.1 

5.0 

22.5 

204 

13.7 

10 

11.25 

n/a 

n/a 

n/a 

7.1 

5.8 

3.7 

11 

46.00 

31.6 

27.3 

9.3 

45.5 

42.5 

25.8 

12 

38.75 

9.2 

5.1 

10.0 

19.7 

13.9 

15.9 

13 

64.25 

8.6 

6.7 

4.9 

15.2 

13.  8 

15.6 

14 

79.50 

n/a 

n/a 

n/a 

37.5 

27.0 

73.9 

15 

35.50 

17.2 

14.5 

6.3 

30.4 

25.4 

20.9 

16 

64.25 

15.2 

11.3 

5.0 

25.2 

13.6 

17.3 

17 

72.50 

13.8 

13.0 

4.2 

26.9 

28.9 

18.1 

18 

68.00 

13.7 

0.3 

1.2 

30.4 

1.3 

4.2 

19 

6.75 

n/a 

n/a 

n/a 

4.7 

2.0 

1.2 

20 

46.00 

15.4 

10.8 

1.2 

20.4 

12.8 

7.7 

21 

38.75 

14.7 

11.4 

1.9 

33.6 

38.3 

19.8 

23 

71.25 

15.9 

2.3 

1.7 

35.1 

5.7 

5.5 

24 

37.25 

19.8 

21.4 

7.5 

40.1 

38.4 

10.5 

25 

38.75 

14.6 

17.9 

16.3 

23.3 

23.3 

19.5 

26 

61.00 

16.8 

15.1 

1.5 

34.8 

37.7 

31.8 

27 

5.25 

n/a 

n/a 

n/a 

7.7 

0.6 

7.7 

28 

26.75 

5.0 

3.0 

3.0 

5.7 

2.9 

3.1 

29 

64.50 

16.4 

13.2 

6.1 

27.4 

14.1 

17.0 

Highlighted  rows  correspond  to  the  sonobuoys  pictured  in  Figs.  16-18. 
Italics  indicate  those  buoys  with  a  lifespan  less  than  24  h. 


the  utility  of  blending  Lagrangian  data  with  model  velocities  and 
not  to  compare  competing  methods  or  determine  their  optimal 
implementations. 

Overall,  the  results  show  that  the  use  of  Lagrangian  data  to  cor¬ 
rect  model  velocity  fields  is  a  powerful  tool  to  improve  trajectory 
prediction.  This  suggests  that  in  many  practical  applications,  such 
as  search  and  rescue  or  pollutant  release  problems,  prediction 
could  be  improved  by  performing  ad  hoc  launches  of  drifters  and 
using  the  data  to  correct  model  results. 

We  stress  that  the  velocity  reconstruction  as  tested  here  has 
limited  forecasting  capabilities.  True  forecasts  with  Lagrangian 
assimilation  require  running  the  model  forward  starting  from  a 
corrected  state,  where  the  mass  variables  are  also  modified,  consis¬ 
tently  with  the  velocity  field.  Early  examples  of  true  forecasts  from 
Lagrangian  assimilations  using  Argo  float  data  in  the  interior  ocean, 
correcting  the  mass  field  based  on  dynamic  principles  in  concert 
with  mass  conservation  (Ozgbkmen  et  al.,  2003),  have  been  re¬ 
ported  in  Taillandier  et  al.  (2006a,  2010).  In  this  context,  the  use 
of  a  dynamically  evolving  correlation  matrix  between  drifters, 
sonobuoys,  and  model  velocities,  as  in  the  Extended  or  Ensemble 
Kalman  filter,  could  provide  further  improvements  in  the  forecast. 
In  addition,  an  approach  that  computes  covariances  dynamically 
may  be  very  useful  for  direct  assimilation  of  sonobuoy  trajectories. 
While  desirable  in  specific  applications,  incorporating  sonobuoy 
trajectories  is  expected  to  be  more  challenging  than  those  of  Argo 
floats  and  drifters,  since  they  sample  a  significant  portion  of  the 
upper  ocean  rather  than  a  single  level. 

Several  related  issues  have  to  be  considered  regarding  assimila¬ 
tion  of  Lagrangian  data  sampling  the  upper  ocean.  First,  Lagrangian 
sensors  might  sample  submesoscale  flows  that  can  be  present  in 
the  surface  layers,  but  that  might  not  be  correctly  resolved  by 
the  operational  models.  This  raises  the  question  of  whether  this 
data  should  be  filtered  to  correct  only  the  mesoscale  component 
(as  done  for  LAVA  in  the  present  paper),  or  submesoscale  signa¬ 


tures  should  be  retained  or  possibly  used  to  improve  model  param- 
eterizations.  A  second  issue  is  that  a  simple  geostrophic  balance  is 
not  expected  to  hold  and  cannot  therefore  be  used  to  correct  the 
mass  field  as  in  Ozgbkmen  et  al.  (2003).  Upper  ocean  and  mixed 
layer  dynamics  imply  the  presence  of  Ekman  dynamics  as  well  as 
possible  strong  nonlinearities  related  to  fronts  and  submesoscale 
features.  It  seems  to  us  that  more  knowledge  is  needed  of  such 
dynamics  to  guide  a  meaningful  state  correction  for  sequential 
data  assimilation.  Alternatively,  this  dynamic  information  could 
be  incorporated  implicitly  in  a  forecasting  system  that  computes 
the  correlations  as  they  evolve,  for  example  from  an  ensemble 
analysis.  Such  an  approach  appears  promising  in  scenarios  where 
the  dynamics  are  nonlinear  and  highly  complex. 
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Appendix  A.  Lagrangian  variational  analysis  (LAVA)  method 

Drifter  data  positions  tote  and  model  velocity  fields  Umod  are 
blended  together  to  reconstruct  two-dimensional  flows.  The  pur¬ 
pose  of  the  methodology  is  to  correct  circulation  structures  whose 
size  and  dynamics  can  be  well  represented  by  the  model  but  may 
appear  shifted  in  space  and  time  relative  to  observations  due  to 
inaccuracies  in  model  initialization  and/or  forcing.  The  circulation 
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structures  to  be  corrected  are  characterized  a  priori  to  have  a 
length  scale  R,  a  time  scale  T,  and  a  Lagrangian  time  scale  Tj,. 

Considering  such  scaling  {R,T,Ti)  constant  over  a  well-defined 
domain  (shelf  or  offshore),  the  reconstruction  consists  of  comput¬ 
ing  time-independent  sequential  velocity  corrections  Au(to)  for 
successive  sampling  sequences  [to,to  +  t]  in  the  neighborhood  of 
drifter  trajectories.  The  sequence  duration  t  is  chosen  to  be 
approximately  Tj.  for  a  consistent  sampling  of  the  corrected  struc¬ 
tures.  The  time-dependent  velocity  field  is  then  reconstructed  as 

Ue(t)  =Um„d(f)  +  Au(to)  for  te  [to-T/2,to  +  T/2].  (A.l) 

For  each  sequence  [to, to  +  t],  an  inverse  problem  is  solved.  The 
velocity  correction  Au  minimizes  the  misfit  between  trajectories 
simulated  using  Umod  and  obsei'ved  positions  Tots  (Taillandier 
et  al.,  2006a).  This  minimization  problem  Is  solved  using  a  varia¬ 
tional  approach  formulated  as  follows: 


Calculation  of  the  vertical  component  of  relative  vorticity  from 
Eq.  (B.l)  in  a  Cartesian  coordinate  system  leads  to  a  Flelmholtz 
equation  for  W.  Here,  '¥  is  represented  as  a  set  of  eigenfunctions, 
which  we  call  Dirichlet  modes  (i/'n)-  The  ij/„,  which  may  be  thought 
of  as  streamfunction  or  vorticity  modes  with  zero  horizontal  diver¬ 
gence,  are  solutions  to 

VVn  +  ^lAn  =  0,  lAnlao  =  0.  (B-2) 

where  D  denotes  the  domain  and  dD  its  boundary. 

Calculation  of  the  vertical  velocity  component  from  Eq.  (B.l)  in 
a  Cartesian  coordinate  system  leads  to  a  Helmholtz  equation  for 
<P  =  d^jdz.  We  represent  cP  as  a  set  of  eigenfunctions,  which  we 
call  Neumann  modes  {(pm)-  The  ipm,  which  may  be  thought  of  as 
velocity  potential  or  divergence  modes  with  zero  relative  vorticity, 
are  solutions  to 


Stage  1:  The  position  robs(  to  "t)  is  predicted  from  a  trajectory 
simulation  solving  the  nonlinear  differential  equation 

dtr  =  Umod(r(t),t)  for  t6[to,to  +  T],  withr(to)  =  rote(to), 

(A.2) 

where  d[  is  the  first  order  derivative  in  time.  This  position 
prediction  can  be  written  as  r(to  +  t)  =  Hivt(Umod)  for  the 
nonlinear  operator  H^l-  The  distance  between  observed 
and  simulated  positions  is  then  expressed  by  the  cost 
function 

1  T 

J  =  2  (rote(to  +  t) - HNi(Umod))  (FoOsjto  +  T)  -  HNi(Umod)), 

(A.3) 

where  ^  denotes  the  vector  transpose.  Components  of  the 
position  difference  are  assumed  to  be  independent.  They 
are  associated  with  Gaussian  homogeneous  errors,  so  that 
the  observational  error  covariance  matrix,  used  for  the 
scalar  product  in  Eq.  (A.3),  is  the  identity  matrix. 

Stage  2:  The  optimal  velocity  correction  Au  is  estimated  by  min¬ 
imizing  J  using  a  steepest  descent  procedure,  along  the 
gradient 

VJ  =  -BH’'(roj,s(to  +  T)  -  HNi(Umod)),  (A.4) 

where  B  is  the  background  error  covariance  matrix,  H  the 
tangent  linear  operator  associated  with  HNL(Umod).  and 
its  adjoint  operator.  B  is  built  by  finite  iterations  of  the  dif¬ 
fusion  equation  (Derber  and  Rosati,  1989;  Weaver  and 
Courtier,  2001)  in  order  to  spread  each  along-trajectory 
velocity  correction  throughout  its  neighborhood.  Note  that 
the  number  of  iterations  depends  on  the  characteristic 
length  scale  R.  H  is  defined  by  the  perturbation  equation 
associated  with  Eq.  (A.2),  as  explicitly  derived  in  a  discrete 
formalism  by  Taillandier  et  al.  (2006a). 


(B.3) 


where  n  is  the  unit  outward  normal  vector  on  the  boundary. 

From  Eq.  (B.l)  and  restricting  out  attention  to  the  horizontal 
velocities,  the  gradients  of  ip„  and  (jjn,  can  be  expressed  as 


-QlAn  dljjA 
dy  '  dx)' 


(B.4) 


d4>m\ 
\dx  '  dy  )' 


(B.5) 


The  analysis  domain  used  here  is  rectangular,  so  that  the  \j/n 
and  ipm  and  their  gradients  are  known  analytically.  Both  eigen¬ 
function  sets  are  ordered  by  decreasing  spatial  scale,  with  the 
first  mode  having  the  largest  spatial  scale.  Note  that  both  the 
ipn  and  (fm  result  in  no  normal  flow  on  the  domain  boundary. 
Since  our  analysis  domain  is  in  the  open  ocean,  flow  across  the 
domain  boundary  is  anticipated.  To  account  for  the  normal  com¬ 
ponent  of  the  flow  at  the  domain’s  open  boundaries,  a  boundary 
velocity  potential  solution  0  Is  calculated  numerically  (using  a 
Matlab  implementation  of  the  generalized  minimum  residual 
method)  as  the  solution  to 


V  0  —  Sg,  (n  •  V0)  Ifljj  —  (n  •  Umod)  laD> 


(B.6) 


where  Umod  are  model  velocities  and  Sq  is  a  source  term  that  ac¬ 
counts  for  the  net  flow  into  the  domain  through  its  open  bound¬ 
aries,  defined  as 


c  foD^'  ^moddl 
ff.dxdy  ■ 

From  0,  boundary  solution  velocities  are  computed  as 
Utdry  =  V0(X.y). 


(B.7) 


(B.8) 


Appendix  B.  Normal  mode  analysis  (NMA)  method 

The  description  here  follows  Llpphardt  et  al.  (2000).  A  more  de¬ 
tailed  derivation  can  be  found  in  Eremeev  et  al.  (1992a,b)  and  Chu 
et  al.  (2003). 

The  three-dimensional  incompressible  velocity  field  Is  ex¬ 
pressed  in  terms  of  two  scalar  potentials  as 

u  =  Vx  [!{(-?')  + Vx  (B.l) 

where  k  is  the  unit  vector  in  the  vertical  direction.  Note  that  this 
form  ensures  that  the  velocity  field  is  exactly  incompressible  In 
three  dimensions. 


For  convenience,  we  will  refer  to  a  single  set  of  basis  functions 
Uii  that  includes  both  the  u°  and  uj(.  With  this  notation,  the  com¬ 
plete  estimated  velocity  field  at  a  given  depth  Is  given  by 

K 

Ue(x,y,t)  =  UMo,(x,y,t)  +  ^a,i(t)Ufc(x,y).  (B.9) 

k=l 

At  each  time  t,  the  amplitudes  a);(t)  are  solutions  to  the  minimiza¬ 
tion  problem 

mm||Aa-b||2,  (B.IO) 

where 
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a2(t) 

? 

.Mt). 

U,(X,)  Uk(X,) 

U,(X2)  •••  U|f(X2) 

.  .  ’ 

••• 

■Umod(Xi,t) -UMo-(Xi,t) 

Umod(X2,t)  -UMo,(X2,t) 
.Umod(XNgri<|i  0  ~  U[i^;j,(Xn^.^,  t)  _ 


(B.ll) 


The  system  is  overdetermined  for  Ngnd  >  K,  and  the  Ok  are  deter¬ 
mined  through  least  squares  minimization.  If  Ngnd  =  K,  then  Aa  =  b 
can  be  solved,  reproducing  the  model  velocities  exactly  on  the  grid. 
When  Ngrid  <  K  the  underdetermined  system  may  be  solved  in  a 
minimum  norm  sense.  The  last  two  cases  do  not  occur  in  this 
study. 

When  drifter  velocities  are  also  available,  they  might  simply  be 
added  as  additional  elements  in  b  above.  For  the  NMA  mappings 
used  here,  however,  there  are  typically  thousands  of  model  veloc¬ 
ities  and  no  more  than  30  drifter  velocities.  In  this  case,  treating 
drifter  velocities  the  same  as  model  velocities  gives  the  drifter 
observations  very  little  weight. 

Since  we  anticipate  differences  between  model  and  drifter 
velocities  and  we  desire  to  use  the  drifters  to  correct  the  model 
forecasts,  the  drifter  velocities  are  used  as  a  separate  constraint 
on  the  least-squares  minimization  problem  (B.IO)  (Toner  et  al., 
2001).  Drifter  paths  rp(t)  (p  =  1,. .  ..Ndri/t)  and  the  associated  veloc¬ 
ities  /ip(t)  (p  =  1,. .  ..Ndrift),  computed  from  Hermite  cubic  interpola¬ 
tions  of  the  positions  along  these  paths,  provide  this  additional 
constraint: 

Ba  =  d,  (B.12) 

where 


d  = 


u,(ri(t))  UK(r,(t)) 

u,(r2(t))  UK(r2(t)) 

i«2(r2(t))-UMo,(r2(t)) 


(t))  -UMo-(rN,„,,(f)) 


(B.13) 


In  the  case  at  hand,  the  domain  is  a  rectangle,  leading  to  analytic  ba¬ 
sis  functions  Uk,  which  are  evaluated  directly  to  generate  velocities 
at  drifter  locations  U(((rp(t)).  The  numerical  boundary  solution  Ut,dry 
is  bilinearly  interpolated  to  the  drifter  locations. 

The  least  squares  problem  (B.IO)  subject  to  the  constraint  (B.12) 
has  a  unique  solution  when  B  has  full  row  rank  and 


(B.14) 


has  full  column  rank,  as  is  always  the  case  here.  We  solve  this  system 
using  a  Matlab  implementation  of  the  QR  factorization  technique. 


Appendix  C.  Sonobuoy  field  drift  model  (SFDM) 


SFDM  is  a  collection  of  Matlab  routines  and  a  Fortran  module 
that  calculate  sonobuoy  trajectories  from  their  initial  positions 
and  hydromechanical  properties,  using  a  model  velocity  field. 
Based  on  the  initial  conditions,  a  velocity  vertical  profile  is  devel¬ 
oped  for  each  buoy  by  three-dimensional  linear  interpolation  of 
the  u  and  v  model  velocity  components.  These  profiles  are  passed 
to  the  sonobuoy  drift  response  model  to  derive  the  drift  velocity, 
which  is  in  turn  used  to  update  the  sonobuoy  position  for  the  next 
time  step  with  a  simple  forward  Euler  step.  The  code  can  be  found 
in  Hammond  (2005). 

The  drift  response  model  is  a  modified  version  of  the  Navy  stan¬ 
dard  sonobuoy  model  FF2E,  described  by  Wang  and  Moran  (1980). 
FF2E  is  a  two-dimensional  steady  state  cable  model  that  is  used  as 
a  standard  design  and  evaluation  tool  in  the  commercial  sonobuoy 
industry.  It  predicts  the  steady  state  response  (including  drift 
velocity)  of  a  free-floating  cable-body  system  to  a  two-dimensional 
current  profile.  It  has  been  extensively  tested  and  validated  (e.g., 
McEachern,  1980).  Wind  drag  is  incorporated  into  EF2E,  but  was 
not  used  in  the  present  application  (see  Section  4). 

Real  world  current  fields  typically  have  a  complex  three-dimen¬ 
sional  structure.  FE2E,  however,  is  restricted  to  two  dimensions,  and 
the  sonobuoy  response  is  calculated  separately  for  two  orthogonal 
horizontal  velocity  components.  To  minimize  the  resulting  drift  er¬ 
ror  of  up  to  nearly  30%  (Hammond,  2005),  these  components  do 
not  necessarily  coincide  with  u  and  but  rather  are  rotated  for  one 
of  them  to  align  with  the  predominant  current  direction,  determined 
by  a  weighted  mean.  The  rotation  angle  a  is  given  by 


a  =  ^  arctan 

Z 


(v(Z)\  U^(z) +  11^(2) 


(C.l) 


EE2E  considers  a  free-floating  system  consisting  of  a  surface 
buoy,  an  arbitrary  number  of  cable  segments  and  intermediate 
bodies,  and  a  terminal  weight  at  the  bottom.  The  ingested  velocity 
profile  is  assumed  to  be  a  function  of  depth  alone.  At  equilibrium, 
the  entire  system  will  drift  at  a  velocity  Up  and  the  surface  huoy 
will  have  a  draft  depth  H.  These  two  variables  are  determined  iter¬ 
atively.  Eor  each  iteration,  the  current  values  for  Uq  and  H  are  used 
to  derive  the  tension  and  angle  on  the  cable  just  below  the  buoy. 
The  equilibrium  equations  are  then  integrated  down  the  cable, 
accounting  for  the  intermediate  bodies  along  the  way.  Equilibrium 
is  checked  at  the  bottom  unit: 


eH  =  T„  -  Db,  (C.2) 

ev  =  Tv-  Wb,  (C.3) 

where  and  ey  are  the  horizontal  and  vertical  imbalances,  Th  and 
Ty  are  the  horizontal  and  vertical  components  of  the  tension  just 
above  the  bottom  unit,  and  Db  and  Wb  are  the  drag  and  weight  in 
water  of  the  bottom  unit.  Various  criteria  on  the  imbalances  can 
be  imposed;  here  the  iteration  is  stopped  when  both  of  the  follow¬ 
ing  inequalities  are  satisfied: 

\Sh\  <  ^i\Db\, 

|ev|  <  e2|WB| 

for  ei,2  small  positive  constants  (taken  to  be  0.02  here). 

To  start  the  iteration,  H  is  a  function  of  the  buoyancy  required  to 
support  the  weight  of  the  entire  system  in  water  in  the  absence  of 
currents,  while  Uo  is  determined  by  a  bisection  scheme  to  reduce 
Em  alone.  If  at  the  end  of  this  step  ey  is  still  too  large,  a  simultaneous 
scheme  is  used.  If  convergence  fails  under  this  scheme,  a  solution  is 
found  with  a  staggered  scheme. 


(C.4) 

(C.5) 
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For  the  simultaneous  scheme,  and  H  are  updated  using 


Aud 


AH  = 


5SH\eH\eH 
EMu 
-S\ev\ev 
EM,  ’ 


and 


(C.6) 

(C.7) 


where  E  =  is  the  combined  equilibrium  error,  M,  and  M, 

are  numerical  factors  defined  below,  and  S  and  Sh  are  positive  pro¬ 
portionality  quantities  chosen  to  facilitate  convergence.  The  numer¬ 
ical  factors  are  given  by 


Mu  =  pCdAt^, 

(C.8) 

M,^pgn^, 

(C.9) 

with  p  the  fluid  density,  Ac  the  difference  between  maximum  and 
minimum  current,  g  the  gravitational  constant,  ds  the  diameter  of 
the  surface  buoy,  and  CdAj  defined  by 


CdAt  =  I  CdAs  +  ^2  ^Dkiakdok  +  ^  CdAi 

)<=i 


(CIO) 


For  the  above  equation  the  following  definitions  were  used:  CqAs  is 
the  drag  area  of  the  surface  buoy.  K  is  the  total  number  of  cable  seg¬ 
ments.  Cok,  Lk,  dok  are  the  drag  coefficient,  reference  length,  and  ref¬ 
erence  diameter  of  the  kth  cable  segment,  respectively.  is  the 
drag  area  of  the  kth  body. 

The  two  proportionality  quantities  are  initialized  as  1.  If  Bh 
changes  too  fast  (slow),  Sh  is  increased  (decreased).  If  the  updated 
values  for  H  and  Uo  lead  to  an  increased  error  £,  Aud  and  AH  are 
recalculated  with  the  previous  values  and  a  reduces  value  for  <5.  S 
is  continually  reduced  until  one  of  the  following  three  conditions 
holds  (where  updated  quantities  are  represented  with  a  prime): 


1.  £'<£, 

2.  e'f^/en  >  I  and  \e'„\  >  |e(,|,  where  /  is  the  drag  force  per  unit  length 
normal  to  the  cable, 

3.  5  is  less  than  a  prescribed  lower  bound  (0.2  for  the  first  ten  iter¬ 
ations,  0.05  thereafter). 


In  any  of  these  cases,  <5  is  reset  to  1. 

If  the  simultaneous  scheme  as  described  above  does  not  con¬ 
verge,  a  staggered  scheme  is  used.  This  scheme  searches  for  appro¬ 
priate  values  for  Ud  and  H  successively.  Uq  is  updated  using  a 
bisection  method,  with  the  sign  of  Auq  determined  by  Eq.  (C.6). 
H  is  updated  using 

AH  =  ^.  (C.11) 

The  proportionality  quantity  Sy  is  initially  set  to  0.6  and  adjusted  if 
convergence  is  too  slow  or  changes  in  H  are  too  large. 

Details  of  the  drag  calculations  can  be  found  in  McEachern 
(1980). 
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